// FE_solver_2D_Parabolic.cpp : This file contains the 'main' function. Program execution begins and ends there.
//

#include <iostream>
#include "FEM.h"
#include "TimeDiscrete.h"
int ex_index = 1;
//coefficients
double coe_function(CPoint p,double t)
{
	return 2.0;
}

double load_function(CPoint p,double t)
{
	if (1 == ex_index)
	{
		//example 1 in chapter 4
		return  -3*exp(p.getX()+p.getY()+t);
	}
	else if (2 == ex_index || 3 == ex_index)
	{
		//example 2 in chapter 3
		return -2.0*exp(p.getX() + p.getY());
	}
	return 0;

}
double left_bound_value(CPoint p,double t)
{
	if (1 == ex_index)
	{
		//example 1 in chapter 4
		return exp(p.getY()+t);
	}

	else if (2 == ex_index || 3 == ex_index)
	{
		//example 2 in chapter 3
		return exp(-1 + p.getY());
	}
	return 0;

}
double right_bound_value(CPoint p,double t)
{
	if (1 == ex_index)
	{
		//example 1 in chapter 4
		return exp(2+p.getY()+t);
	}

	else if (2 == ex_index || 3 == ex_index)
	{
		//example 2 in chapter 3
		return exp(1 + p.getY());
	}
	return 0;

}
double bottom_bound_value(CPoint p,double t)
{
	if (1 == ex_index)
	{
		//example 1 in chapter 4
		return exp(p.getX()+t);
	}

	else if (2 == ex_index)
	{
		//example 2 in chapter 3
		return -exp(p.getX() - 1);
	}
	else if (3 == ex_index)
	{
		return 0;
	}
	return 0;

}
double up_bound_value(CPoint p,double t)
{
	if (1 == ex_index)
	{
		//example 1 in chapter 4
		return exp(p.getX()+1+t);
	}

	else if (2 == ex_index || 3 == ex_index)
	{
		//example 2 in chapter 3
		return exp(p.getX() + 1);
	}
	return 0;

}

//nalba_u+ru=q
double r(CPoint x)
{
	return 1;
}


double analytical_solution(CPoint p, int der_x, int der_y,double t)
{
	
	double x = p.getX();
	double y = p.getY();

	if (ex_index == 1)
	{
		//example 1 in chapter 4
		if (der_x == 0 && der_y == 0)
		{
			return exp(p.getX()+p.getY()+t);

		}
		else if (der_x == 1 && der_y == 0)
		{
			return exp(p.getX() + p.getY() + t);
		}
		else if (der_x == 0 && der_y == 1)
		{
			return exp(p.getX() + p.getY() + t);
		}
		else if (der_x == 1 && der_y == 1)
		{
			return exp(p.getX() + p.getY() + t);
		}
	}

	else if (ex_index == 2 || ex_index == 3)
	{
		//example 2 and 3 in chapter 3
		if (der_x == 0 && der_y == 0)
		{
			return exp(x + y);

		}
		else if (der_x == 1 && der_y == 0)
		{
			return exp(x + y);
		}
		else if (der_x == 0 && der_y == 1)
		{
			return exp(x + y);
		}
		else if (der_x == 1 && der_y == 1)
		{
			return exp(x + y);
		}
		else
			return 0;
	}
	return 0;


}
double func_one(CPoint p, double t)
{
	return 1.0;
}
double initial_value(CPoint p)
{
	return exp(p.getX() + p.getY());
}
int main()
{
	
	CFEM fem;
	fem.Init();

	//read mesh from file
	//the mesh file contains the P,T,Pb,Tb,boundary nodes,boundary edges, and the type of the element
	fem.readMesh("..\\mesh\\cha4_exa1_qua.neu");

	//the basis type should be consistent with the mesh file
	BASIS_TYPE basis_trial_type = TWO_DIM_QUADRATIC;
	BASIS_TYPE basis_test_type = basis_trial_type;

	//assign boundary nodes value by the boundary index defined in the .neu mesh file
	//so the assignment of the boundary must be consistent with the mesh file
	//example 1 in chapter 4
	//boundary values are all defined by functions, it depends on the coordinates of the nodes
	fem.assignBoundNodesValue(0, Dirichlet, left_bound_value);
	fem.assignBoundNodesValue(1, Dirichlet, bottom_bound_value);
	fem.assignBoundNodesValue(2, Dirichlet, right_bound_value);
	fem.assignBoundNodesValue(3, Dirichlet, up_bound_value);

	//time discrete
	double time_left = 0.0, time_right = 1.0;
	int number_of_time_step = 10;
	double theta = 0.5;
	CTimeDiscrete time(time_left, time_right, number_of_time_step,theta);

	//assemble mass matrix
	SpMat M = fem.assembleMatrix2D(func_one, 0, basis_trial_type, 0, 0, basis_test_type, 0, 0);
	
	//assemble stiffness matrix	
	//call assemble function for x direction
	SpMat A1 = fem.assembleMatrix2D(coe_function, 0,basis_trial_type,1, 0, basis_test_type,1, 0);//x

	//call assemble function for y direction
	SpMat A2 = fem.assembleMatrix2D(coe_function, 0,basis_trial_type,0, 1, basis_test_type, 0, 1);//y
	SpMat A = A1 + A2;

	Rsv x0=fem.generate_Initial(initial_value);
	Rsv x;
	
	for (int m=0;m<time.getTimeStepNumber()-1;++m)
	{
		double dt = time.getTimeStep(m);
		double tm = time.getTime(m);
		double tmp1 = time.getTime(m + 1);
		Rsv bmp1 = fem.assembleVector2D(load_function, tmp1, basis_test_type, 0, 0);
		Rsv bm = fem.assembleVector2D(load_function, tm, basis_test_type, 0, 0);

		//dt may not be the same for all time step
		SpMat Atilde = M / dt + time.getTheta()*A;
		SpMat Afixed = M / dt - (1 - time.getTheta())*A;

		Rsv btilde = time.getTheta()*bmp1 + (1 - time.getTheta())*bm + Afixed * x0;

		//treat boundary condition
		fem.treatBoundaryConditions(coe_function, tmp1, Atilde, btilde);

		//solve the linear equations
		x = fem.solveLinearEqua(DIREC_PardisoLDLT, Atilde, btilde);

		x0 = x;
	}	
	
	//export result
	fem.Export2VTK("result.vtk", x);

    //error measurement
	fem.errorMesu(2, x, analytical_solution,time.getTime(time.getTimeStepNumber()-1));	

}
